@rhoprof
restore,'grid.sav'

nx=(size(ph))[1]
ny=(size(th))[1]
nz=(size(r))[1]

lev='03'

if(lev eq '04') then begin 
 if (nz eq 47)  then alt=43
 if (nz eq 93)  then alt=90
 if (nz eq 188) then alt=180
 if (nz eq 100) then alt=97
endif
if(lev eq '03') then begin 
 if (nz eq 188) then alt=86
 if (nz eq 100) then alt=86
endif
print,lev,alt

file1='uxy'+lev
file2='vxy'+lev
file3='wxy'+lev

uu = dblarr(nx,ny)
vv = dblarr(nx,ny)
ww = dblarr(nx,ny)

openr,1,file1,/f77
openr,2,file2,/f77
openr,3,file3,/f77
ns = 0L
while (not eof(1)) do begin
  readu,1,uu,t
  if (ns LE 1) then begin
    u = uu 
    print,'Reading files: ',file1
  endif else begin
    u = [[[u]],[[uu]]]
  endelse
  ns ++
end
close,1

ns = 0L
while (not eof(2)) do begin
  readu,2,vv,t
  if (ns LE 1) then begin
    v = vv 
    print,'Reading files: ',file1
  endif else begin
    v = [[[v]],[[vv]]]
  endelse
  ns ++
end
close,2

ns = 0L
while (not eof(3)) do begin
  readu,3,ww,t
  if (ns LE 1) then begin
    w = ww 
    print,'Reading files: ',file1
  endif else begin
    w = [[[w]],[[ww]]]
  endelse
  ns ++
end
close,3

print,size(u)
print,size(v)
print,size(w)

v1=u
v2=v
v3=w

ti=10

 vx = fltarr(nx,ny,ns-1)
 vy = fltarr(nx,ny,ns-1)
 vz = fltarr(nx,ny,ns-1) 
 vtx = fltarr(nx,ny)
 vty = fltarr(nx,ny)
 vtz = fltarr(nx,ny) 
 re_vtx = fltarr(nx,ny)
 re_vty = fltarr(nx,ny)
 re_vtz = fltarr(nx,ny) 
 im_vtx = fltarr(nx,ny)
 im_vty = fltarr(nx,ny)
 im_vtz = fltarr(nx,ny) 
 sp = fltarr(nx/2,nx/2)
 mean_sp = fltarr(nx/2,ns-1)

m=0L
Lx=max(ph*!pi/180.)
Ly=max(th*!pi/180.)
kx=shift(-(nx)/2+indgen(nx),+(nx+1)/2)
ky=shift(-(ny)/2+indgen(ny),+(ny+1)/2)
;nk=fix(sqrt(((nx)/lx)^2+((ny)/ly)^2 )*lx/2)
nk=fix(sqrt(((nx))^2)/2)
;nk=nx/2
spectrum=fltarr(nk,ny,ns-1)
ki=intarr(nk)

for in=0,ns-2 do begin
; vx[*,*,in] = sin(ph[*]*!pi/180.)#cos(4.*th[*]*!pi/180.); v1[*,*,in]
 vx[*,*,in] = v1[*,*,in]
 vy[*,*,in] = v2[*,*,in]
 vz[*,*,in] = v3[*,*,in]

 for iky=0,ny-1 do begin
  vtx[*,iky] = fft(vx[*,iky,in],/double)
  vty[*,iky] = fft(vy[*,iky,in],/double)
  vtz[*,iky] = fft(vz[*,iky,in],/double)

  re_vtx[*,iky]=real_part(vtx[*,iky])
  re_vty[*,iky]=real_part(vty[*,iky])
  re_vtz[*,iky]=real_part(vtz[*,iky])

  im_vtx[*,iky]=imaginary(vtx[*,iky])
  im_vty[*,iky]=imaginary(vty[*,iky])
  im_vtz[*,iky]=imaginary(vtz[*,iky])
 endfor

it=0D
for iky=0,ny-1 do begin 
;iky=ny/2
  for ikx=0,nx-1 do begin
   k=kx[ikx]
   if(k gt 0 AND k le nk-1) then begin
    ki[k]=k
    prods=0.5*((re_vtx[ikx,iky]^2+im_vtx[ikx,iky]^2)+$
	       (re_vty[ikx,iky]^2+im_vty[ikx,iky]^2)+$
 	       (re_vtz[ikx,iky]^2+im_vtz[ikx,iky]^2))
    spectrum[k,iky,in]=spectrum[k,in]+prods
   endif
   it++	
  endfor
endfor


endfor

psp=total(spectrum,3)/ns
psp=total(psp,2)/(ny/2.)
;psp=rh[alt]*total(mean_sp,2)/(ns-1)
;psp=total(spectrum,2)/ns

save,filename='power'+lev+'.sav',psp,ki
plot,ki,psp,/xlog,/ylog,xr=[1,400],xtitle="m",ytitle="E!dK!n(k)",xstyle=1

; line with slope -5/3 
y = 0.95e4*ki^(-5./3.)
 
oplot,ki,y,linestyle=2

end
